﻿#include "DILUPrecond.hpp"

template <typename Cmpt>
void UNAP::DILUPrecond<Cmpt>::calcReciprocalD(Vector<Cmpt> &rD) {
  Cmpt *rDPtr = rD.data();

  forAll(i, nFaces_) {
    rDPtr[uPtr_[i]] -= upperPtr_[i] * lowerPtr_[i] / rDPtr[lPtr_[i]];
  }

  reciprocal(rD);
}

template <typename Cmpt>
void UNAP::DILUPrecond<Cmpt>::precondition(Vector<Cmpt> &w,
                                           const Vector<Cmpt> &r) const {
  Cmpt *wPtr = w.data();
  const Cmpt *rDPtr = rD_.data();
  const label nFacesM1 = nFaces_ - 1;
  label sface;

  pointwiseMul(w, rD_, r);

  forAll(i, nFaces_) {
    sface = losortPtr_[i];
    wPtr[uPtr_[sface]] -=

        rDPtr[uPtr_[sface]] * lowerPtr_[sface] * wPtr[lPtr_[sface]];
  }

  for (label face = nFacesM1; face >= 0; face--) {
    wPtr[lPtr_[face]] -=
        rDPtr[lPtr_[face]] * upperPtr_[face] * wPtr[uPtr_[face]];
  }
}

template class UNAP::DILUPrecond<scalar>;
